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ABSTRACT 



We describe an approach for finding the eigenfrequencies of solar acoustic modes 
jy^ • (p modes) in a convective envelope in the WKB limit. This approximation restricts 

us to examining the effects of fluid motions which are large compared to the mode 
wavelength, but allows us to treat the three-dimensional mode as a localized ray. The 
method of adiabatic switching is then used to investigate the frequency shifts resulting 
from simple perturbations to a polytropic model of the convection zone as well as from 
^r) [ two basic models of a convective cell. We find that although solely depth-dependent 

perturbations can give frequency shifts which are first order in the strength of the 
perturbation, models of convective cells generate downward frequency shifts which are 
second order in the perturbation strength. These results may have implications for 



resolving the differences between eigenfrequencies derived from solar models and those 
found from helioseismic observations. 



6 

^ 1. Introduction 



Comparisons between eigenfrequencies of solar acoustic modes (also known as p modes, since 
pressure provides the restoring force) determined from solar models and eigenfrequencies measured 



via helioseismology show that our understanding of the Sun's structure is incomplete QGough e\ 



al. 1996). In particular, the model frequencies show the largest discrepancies in three distinct 
regions: near the solar core, at the tachocline (the boundary between the radiative interior and 
the convective envelope), and near the solar surface. In this work, we focus on the discrepancy 
near the solar surface, which is thought to be partially a consequence of the incomplete modeling 
of convective effects in solar models. 

In standard solar models, the stratification of the convection zone is determined by 
mixing-length theories in which one lengthscale is used to parameterize all convective motions. 
In reality, the outer portions of the convection zone harbor turbulent convection on scales at 
least as small as granules (1 Mm) and as big as supergranules (20-30 Mm). In addition, cells 
spanning the entire convection zone (200 Mm) may also exist. These turbulent motions affect 
the frequencies and linewidths of p modes, but in a manner which is not yet fully understood. 
Most authors (for example, Brown 1984, Lavely & Ritzwoller 1993, Gruzinov 1998| ) treat the 
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lengthscales of convective motions as small compared to the mode wavelengths. By parameterizing 
the effects of convection as a turbulent pressure, they are able to average over the motions of 
individual cells, thereby simplifying the analysis. Our work lies at the opposite limit. We assume 
that the wavelengths of the modes are small compared to the lengthscale of the convection. This 
formulation allows us to work in the WKB limit and treat the modes as rays following trajectories 
identical to those of classical particles. 

Working in the WKB limit also allows us to use the formalism of Hamiltonian mechanics 
to find the p-mode eigenfrequencies given a model of the convective envelope. We introduce 
two related methods of doing so, EBK (Einstein-Brillouin-Keller) quantization and adiabatic 
switching. After showing that adiabatic switching is better suited for finding the eigenfrequencies 
of complex systems, we apply it to a few simple examples of convective envelopes. 

We find that for convective envelopes the fractional frequency shift scales as the square of the 
ratio of the convective velocity v to the sound speed c. The frequency shift averages to zero to first 
order in v/c and is downshifted to second order because the ray spends more time in the region 
of negative Doppler shift ( Brown 1984 ). A downshift is produced by temperature fluctuations for 



similar reasons. Roughly, 

- 5uj/lu ~ v 2 /c 2 (1) 

where 6u is the change in the eigenfrequency u. Evaluating c at the midpoint of a p mode cavity 
gives 

v 2 

-Suj/uj ^~(wi/400) 2 ^ (2) 

gR& 

where v\ is the convective velocity in km/s and I is the mode degree. For typical velocities of 
0.3 km/s at a depth of 10 Mm (the cavity midpoint of a n = 1, I = 200 p mode), the fractional 
frequency shift is ~ 10~ 4 . Larger £ values will produce larger shifts. Observational determinations 
of p-mode eigenfrequencies quote error bars of 5lo/uj ~ 5 x 10 -4 , suggesting that this mechanism 
could produce a detectable shift. 

We develop the problem, present the theoretical justification for the ray approximation, and 
provide necessary concepts from Hamiltonian theory in §|2[ In §^ we describe two methods for 
finding the eigenfrequencies of convective envelopes, EBK quantization and adiabatic switching. 
Our results for three convective structures, a vertical sound speed perturbation and two simple 
models of convective cells, are presented in §||. In §|5| we discuss our results and outline several 
possible extensions of the work. 



2. Theoretical Background 



Like an organ pipe, the Sun has resonant modes of oscillation (see Gough &: Toomre 1991 
for a review). Each eigenmode is denoted, in a manner similar to the designation of orbitals in 
the hydrogen atom, by its radial order n, spherical harmonic degree £, and azimuthal order m. 
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Although in general the eigenfrequency uj associated with each mode depends on n, £, and m, the 
assumption of spherical symmetry in a model creates a degeneracy so that u is independent of m. 
In the Sun, rotation breaks the spherical symmetry and hence the degeneracy. 

Solar acoustic modes are radially confined within a resonant cavity by the sound speed 
structure and stratification of the solar interior (the exception being modes with I = which 
penetrate to r = 0). As £ increases at fixed n, the inner turning point of a mode moves outward in 
radius. Although there is some dependence on n, the resonant cavities of p modes with £ ;> 60 lie 
completely within the convection zone and we expect that these modes are most strongly affected 
by convective structures. Since we are primarily interested in the turbulent convection near the 
solar surface, we will consider oscillations with I ^ 200 which are trapped within the outer 10% 
(by radius) of the Sun. 

Instead of treating the complete problem of the effects of convective structures of arbitrary 
scale on the structure and frequencies of high-£ modes, we limit the problem still further. 
Following the lead of time-distance helioseismology ( Duvall et al. 1993| ), we approximate the 



three-dimensional mode as a localized ray, in the process losing any information concerning 
perturbations in the structure of the mode, but allowing the ray to be described as a Hamiltonian 
system. In order to make this simplification — equivalent to the limits of geometrical acoustics 



( Landau &: Lifshitz 1959 ) and the WKB approximation — we must restrict ourselves to considering 



the effects of convective structures which are large compared to the wavelength of the mode. 



2.1. The Modal Equation 

We begin with the linearized equations of hydrodynamics for the vector displacement 5r of a 
fluid parcel. The continuity equation is 

p' + V-(p5r) = 5p + pV-{5r)=0 (3) 

and the momentum equation is 

d 2 5r 

P-Qp = -W + PS' + p'g, (4) 

where df and /' are the Lagrangian and Eulerian perturbations to an equilibrium quantity /, 
and are related through 5f = f + (5r ■ V)/. The density, pressure, and gravity are denoted by 
p, p, and g respectively where boldface symbols denote vector quantities. Throughout this work 
all fluid motions are assumed to be adiabatic, hence the density and pressure perturbations are 
connected by 

6p = —dp, (5) 



P 

where the adiabatic exponent 

'<91np 
dlnp 



(6) 
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with the subscript s indicating that the derivative is taken at constant specific entropy. Finally, 
Poisson's equation relates g' and p': 

V • g' = -4ttG P ', (7) 
where G is the Newtonian gravitational constant. 

Furthermore, in order to simplify the calculations we: 

1. Assume a plane-parallel geometry. We are concerned with high-degree modes which are 
trapped close to the surface and feel the curvature of the Sun only as a small perturbation. 
This assumption allows us to treat the equilibrium gravitational acceleration as a constant: 
g = 

2. Assume the region of propagation takes the form of a two-dimensional slab. The vertical 
coordinate, increasing inwards from the surface, is given by z and the horizontal coordinate 
is x. The horizontal length of the slab, L, is the Sun's circumference: L = 2ttRq. Finally, 
we assume periodic boundary conditions in the horizontal direction so that f\ x =o = f\x=L 
for any variable /. 

3. Neglect perturbations to the gravitational acceleration so that g' = 0. This is known as the 
Cowling approximation and is justified when either n or I are large flCowling 1941 ), the 
latter being the case in this work. 



It has been shown flGough 1993 ) that, given the above assumptions, the system of equations 



©~(0) can De reduced to a single equation for the variable ^ = — p x l 2 b]) = c 2 /? 1//2 V • 5r: 

/a 2 ,\a 2 f 9 5 2 , 9 9, 



where V? is the horizontal Laplacian operator and 



„2 



Tip 



P 

is the adiabatic sound speed. The acoustic cutoff frequency uj c is defined by 

„2 



0) 



^E^ll^zT//,,) (10) 

where 



dz 

is the density scale height. The buoyancy, or Brunt- Vaisala, frequency N is defined by 
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with iV 2 < corresponding to convective instability. 

A complete solution of equation (g) yields both p and g modes, the latter having buoyancy 
as their restoring force. Since the solar convection zone is nearly adiabatically stratified, we have 
iV 2 ~ there. Neglecting the third term in equation (||) so that only p modes are allowed, we 
arrive at a modified wave equation 

B 2 \ 

^2+w 2 U-c 2 V 2 ^ = 0. (13) 
Note that both c and oj c can be functions of position. 



2.2. The Ray Approximation 

At this point, equation (|l3[) may still be solved as an eigenvalue problem for the structure and 
eigenfrequencies of the p modes. In order to consider the problem in the ray limit, we must make 
one further approximation. Consider, for the moment, a simple solar model where oj c = and 
c is a constant. Equation ( |l3| ) then reduces to the Helmholtz equation, which has a plane- wave 
solution 

^ = * e i ( k - r -^). (14) 
The ray associated with this plane wave follows a trajectory perpendicular to the wavefronts. 

Now, consider instead a system where the wavelength of the mode is much shorter than the 
lengthscale over which equilibrium quantities (such as c, p and p) vary. Letting the ratio of these 
lengthscales be denoted by A -1 — which we treat as a constant, small parameter — we write the 
solution to equation ( |l3| ) as a plane wave with varying amplitude, and phase, 

f = * (r,i)e iA *( r >*). (15) 

Although we have not made any approximations thus far, the form of equation (^) suggests that 
^ is rapidly oscillating compared to the background state, implying a solution which may be 
thought of as locally planar. These planar wavefronts allow the motion to be approximated by a 
ray. 

Using this ansatz and following the work of Gough (1993), we expand equation (|l^) and 
equate powers of A. This process is equivalent to the WKB approximation and gives the leading 
equation 



- 



dt ) V A 



c 2 



d<§> 



2 



0. (16) 



Next, making the analogy between equations ( |l4| ) and (|T^), we define the local frequency u> and 
wavenumber k as 

lo = -A— and k = A— (17) 
dt dx v ' 
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These identifications allow us to write equation (16) in the form of a dispersion relation: 

lj = (c 2 k 2 + uj 2 c ) 1 i. (18) 



Furthermore, equation (|l7| ) also implies equations of motion: 

dk dco dx duj 

dt dx dt <9k' 

where 4z denotes the derivative with respect to time along a ray path. 



(19) 



Equations (|l9|) are Hamilton's equations for a Hamiltonian H = to in terms of the canonical 
positions % = X{ and momenta pi = k{. The ray's dispersion relation, cj(x, k, t), is identified as the 
functional form of the Hamiltonian. While this result may appear serendipitous, it is well-known in 



quantum mechanics (Goldstein 1965) that high-frequency solutions to the wave equation (in other 
words, solutions which oscillate rapidly compared to any background variation) follow dynamics 
similar to classical particles. 



2.3. Action variables and adiabatic invariants 

The identification with Hamiltonian mechanics described in the preceding subsection allows 
us to take advantage of the extensive results available in the field. In this subsection, we review 
several results from Hamiltonian theory which will be used later. More details can be found in a 
number of classical mechanics texts (for example, poldstein 1965 , Arnol'd 1978 , Lichtenberg 



Lieberman 1983| , [Tabor 1989|) 



First, we present two results which follow from the structure of the governing equations. In 
any Hamiltonian system it is known that 

dH dH 

Ht = «" < 20 > 

Hence, for systems where the Hamiltonian is not explicitly time-dependent, it is a constant of the 
motion. Second, if the Hamiltonian is cyclic (independent) of one of the canonical positions qt, 
then the corresponding canonical momentum pi is a constant of the motion: 

dH , dpj , . 

— = implies = 0. (21) 

dq { F dt y J 

Certain special Hamiltonians, named integrable systems, have equal numbers of independent 
constants of the motion as degrees of freedom. Such systems allow a change of variables where all 
of the new canonical positions, termed angle variables, are cyclic. The new canonical momenta are 
the actions, 1^, and are defined as 



h 

27T J Ck 



r<f P d ^ ( 22 ) 
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for k = 1, . . . , re, where n is the number of degrees of freedom of the system. The curves Ck are 
topologically-independent closed paths in phase space which do not have to be ray trajectories. 
Since the angles are cyclic, each action is a constant of the motion resulting in re total constants. 
Since a Hamiltonian system with n degrees of freedom has a 2n-dimensional phase space, 
integrability ensures that the phase-space trajectories are confined to (2n — n) = re-dimensional 
surfaces which have the topology of re-dimensional tori. An integrable system is placed on a 
particular invariant torus by its initial conditions and cannot stray from this torus as the system 
evolves in time. 

A concrete example is provided by the Hamiltonian of equation ( |l8|) in a two-dimensional 
geometry. The system has two degrees of freedom, so its phase space has four dimensions. If 
we assume that the Hamiltonian is independent of time, then the frequency uo is a constant. 
Furthermore, if we consider a system where c and oj c are independent of one of the position 
coordinates, say x, then the corresponding wavenumber k x is also a constant of the motion. The 
existence of these two constants ensures the integrability of the system and restricts phase-space 
motions to a two-dimensional surface, topologically equivalent to a two-dimensional torus. 

General Hamiltonian systems, however, are usually non- integrable. The fate of invariant tori 
in non-integrable systems was addressed by the KAM (Kolmogorov-Arnol'd-Moser) theorem. It 
concerns systems which are nearly-integrable in the sense that they may be written as 

H = H (I) + eH 1 {l,®), (23) 

where I and are vectors of the action and angle coordinates, respectively, and e is a small 
parameter. Roughly, the KAM theorem states that for sufficiently small values of e, most invariant 
tori are preserved. In other words, for small perturbations from integrability, most sets of initial 
conditions remain on invariant tori. However, the destroyed tori are distributed throughout the 
invariant ones in phase space and as the perturbation increases, their density increases. For a 
strong enough perturbation, no invariant tori remain. 



3. Methods 

Having shown that within the WKB approximation we may treat p modes as a Hamiltonian 
system, the powerful techniques developed for such systems become available for our use. Instead 
of solving the eigenvalue problem of equation (^) , we seek the eigenfrequencies of the Hamiltonian 
given by equation (|l8|). However, equations (|l^) describe the motion of any ray; we do not know 
a priori which rays correspond to eigenmodes. We consider two methods of isolating the systems 
corresponding to eigenmodes: EBK quantization and the method of adiabatic switching. Both 
of these methods have been used in chemical physics to understand molecular spectra ( Patterson 



1985 , [Shirts et al. 1987 ). Gough (1993) discussed EBK quantization of stellar waves but did not 



fully implement the method. 
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In this paper we are interested in finding the eigenfrequencies of non-integrable systems which 
can be expressed as a perturbation to an integrable system as in equation (p3|). Our base solar 
state, corresponding to Hq in equation ( p3[ ) and which we will show below to be integrable, will be 
a two-dimensional, adiabatically stratified, plane-parallel polytrope. Polytropic systems are those 
where the pressure p and density p are assumed to be related by 



P_ = P_ 

Po \Po 



(24) 



where po and po are constants and p is the polytropic index, which for an adiabatically stratified 
system may be written as p = 1 + 1/T± in terms of the adiabatic exponent defined in equation (|6|). 
The outer portion of the solar convection zone is well-approximated by a polytrope with p ~ 3. 

Enforcing hydrostatic equilibrium leads to the relations 

iUf^r 1 and £.-(±Y (25) 
Po \zoJ Po \zoJ 

where zq is a constant. Finally, the adiabatic sound speed and cutoff frequency are given by 

2 9z , 2 p(p + 2 ) 2 f «<\ 
c = — and lo = - — c (26) 

p \z l 



where g is the gravitational acceleration. 

Figure p] shows a sample raypath in an adiabatically stratified, two-dimensional polytropic 
envelope with periodic horizontal boundary conditions and p = 3. The vertical structure of this 
ray corresponds to a mode with n = 1, £ = 200 (see below for a description of the mode numbers); 
however, the horizontal structure has been altered for legibility - it is equivalent to a mode 
with n = 1,£ = 12. If drawn to scale, the box would be a thin ribbon, with an aspect ratio of 
~ 300 : 1. Furthermore, the ray would fit ~ 80 horizontal wavelengths into the box instead of ~ 5 
as displayed. The ray is confined in the vertical direction by two boundaries, or caustic surfaces, 
where k z = 0. If we continue to trace the ray, it will eventually fill the entire cavity. 



3.1. EBK quantization 



While finding the eigenvalues of a Hamiltonian is usually considered as a quantum-mechanical 
problem, historically it was first approached from the standpoint of classical mechanics. The 
quantization which led to the Bohr model of the hydrogen atom was an early attempt to find 
eigenvalues through the quantization of the action (in this case, the orbital angular momentum). 
Later refinements resulted in EBK semi-classical quantization ( Einstein 1917 , Brillouin 1926| , 
Keller 1958) which states that the action 1^ of equation (p2|) is quantized: 



2^i p " dq= ( nfc + ! r) 



(27) 
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The quantum numbers rik and rrik are non-negative integers, with representing the number of 
caustic surfaces (those surfaces which define the envelope of possible trajectories) crossed by C^. 

As an example of EBK quantization, consider the polytropic envelope described above. 
Christensen-Dalsgaard (1980) showed that the eigenfrequencies for this system are given by 



2 Zgkx ( . V 

to = n H — 

fi V 2 



for 



n 



1,2,... 



(28) 



where k x = £/Rq where Rq is the solar radius and £ is an integer. The latter condition arises from 
requiring that t wavelengths fit around the Sun's circumference. 



EBK quantization can nearly reproduce these results. Combining the expressions for c 2 and 
from equation (26) with the dispersion relation of equation (18), we may solve for k z : 



LO 2 - LO 2 



[ILO 



Az 2 \ ii 



(29) 



However, to is constant since the dispersion relation does not explicitly depend on time. 
Furthermore, k x is also constant since the unperturbed quantities are independent of x. Evaluating 
equation (|27]) on a curve of constant x, and noticing that this curve touches two caustics, one each 
at the top and bottom of the cavity, we have 



I z = - 



7T 



Z2 



k z dz 



1 

n 

2 



for 



ii 



1,2, 



(30) 



where z\ and Z2 are the upper and lower turning points of the ray. With k z given by equation 



|29| ), the integral in equation (30) can be done analytically. The final result is 

2gk x 



LO 



1 a ( 2\ a 



for 



n 



1,2,... 



(31) 



If we also evaluate equation ( J27| ) on a curve of constant z, we see that the curve touches no 
caustics and hence 

1 /•27tRq 

4 = — / k x dx = £ for £ = 0, 1, . . . (32) 



2tt 



This gives: 



k x — 



for 



0,1,... 



(33) 



Equation (|33| ) is the correct horizontal quantization condition. Equations (^g) and (|3l| ) for 
the vertical quantization are not the same, however. Taking n = 1 and [i = 3/2, the relative error 
in the eigenfrequencies is ~ 6%, with the agreement improving as either nor /i increase. The 
magnitude of this error is much larger than the accuracy to which the true p-mode eigenfrequencies 
may be determined. However, we are investigating frequency shifts, not magnitudes, and so the 
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small differences between EBK frequencies and the true frequencies are expected to have only a 
small effect on our results. 

While EBK quantization appears to be a fruitful method for determining the eigenfrequencies, 
it is not generally applicable. Analytically evaluating the quantization conditions of equations 
( |3C)| ) and (^) for all but extremely simple expressions for the sound speed and cutoff frequency is 
difficult. In particular, if these parameters depend on x, the horizontal wavenumber k x is not a 
constant of the motion and the evaluation of the integrals is impractical. 

An alternate approach implements the quantization conditions by propagating the ray 
trajectories and evaluating equations (|30| ) and (^) numerically. By thoughtfully choosing the 
paths Cfc, the EBK equations may be reduced to quantizing the areas of Poincare surface of 
sections (phase space slices in the pi-qi plane). While Noid and Marcus (1975) used this approach 
with some success to find the energy levels of an anharmonically coupled pair of oscillators, it is 
not optimal for many non-integrable systems. As such systems move farther from integrability, 
increasing numbers of invariant tori are destroyed. The initial conditions corresponding to these 
tori are then free to wander in phase space, giving the surfaces of section a characteristic "fuzzy" 
appearance, an early example of which is seen in the work of Henon and Heiles (1964). The form of 
these surfaces of section makes numerical determination of their areas very difficult. Consequently, 
we seek another method for implementing the EBK quantization rules for non-integrable systems. 



3.2. Adiabatic Switching 

Ehrenfest (1917) was among the first to propose the concept of adiabatic switching, although 
his work has been superseded by modern treatments. In essence, his hypothesis was that if 
a system is changed in a reversible, adiabatic way, then allowed motions will be smoothly 
transformed to allowed motions. In the language of § |2.3| , this is analogous to saying that a system 
will remain on an invariant torus if it is altered in a smooth manner. A thorough review of the 
modern implementation of the method of adiabatic switching and its applications can be found in 
Skodje & Cary (1988). 

In adiabatic switching, the system is initialized in an eigenstate of a integrable Hamiltonian 
for which the eigenfrequencies are known. The system is allowed to evolve as the strength of a 
perturbation (not necessarily small, although small in our case) to the original Hamiltonian is 
slowly increased, eventually leaving the system under the direction of a new, usually non-integrable, 
Hamiltonian for which the eigenfrequencies are not known. Due to the adiabatic nature of the 
switch, however, the original eigenstate has slowly relaxed into an eigenstate of the new system 
from which the new eigenfrequency can then be determined. 

While adiabatic switching has been used extensively, for example, to determine the energy 
levels of systems with many degrees of freedom in complicated external electromagnetic fields, it 
does not have a complete mathematical justification. For systems with one degree of freedom, it is 
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known that the action is an adiabatic invariant which remains constant under slow perturbations 
to the Hamiltonian, where "slow" means "long with respect to other characteristic timescales of 



the problem". In the language of £2.3, the initial invariant torus can be deformed into an invariant 
torus of the final state. For systems with more than one degree of freedom, no such theorem 
exists and the overall dynamics are not as well understood. In particular, the KAM theorem 
tells us that some invariant tori of the initial system will be destroyed with the introduction of 
non-integrability. What happens when the system switches through such destroyed tori, found 
everywhere in phase space, is not known in general. 

Despite these considerations, the empirical justification for adiabatic switching is quite 
strong. Skodje & Cary (1988) explore these and other difficulties and conclude that, with some 
precautions, adiabatic switching is an excellent method for implementing EBK quantization in 
non-integrable systems. 

The mathematical formulation of adiabatic switching is fairly straightforward. The non- 
integrable Hamiltonian for which the eigenfrequencies are desired, which in our case is a function 
of position and wavenumber, is written as a sum of two terms 

H(x,k)=H + H 1 , (34) 

where Hq is an integrable Hamiltonian for which the eigenfrequencies are known from EBK 
quantization and Hi is a term which includes everything else. Into equation (|34|), we introduce a 
time-dependent switching function \(t): 

H(x,k) = H + \(t)H 1 (35) 

which satisfies 

A(0) = and A(T) = 1 (36) 

for some time T which, in order to ensure that the transition from Hq to H is adiabatic, is taken to 
be much longer than other timescales associated with the system (for instance, the characteristic 
propagation time t ray rj a; -1 ). In addition, A is chosen such that its first few derivatives are 
continuous at the endpoints. In this work, A is taken to be 

A( t ) = ---sm(-). (37) 

Johnson (1985) has a discussion of the effects of different switching functions. 

At t = our initial conditions are such that we satisfy the EBK quantization conditions for 
the Hamiltonian Hq. We then numerically integrate Hamilton's equations of motion ( |i~9| ) under the 
influence of the Hamiltonian H given in equation (^). As t increases so does A(t), adiabatically 
switching on the non-integrable term H i , and allowing the eigenfrequency of the system to slowly 
adjust from its original value. At t = T we arrive at the eigenfrequency of the full Hamiltonian H. 
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3.3. Numerical Methods 

For our reference state Hq, we use the plane-parallel, adiabatically stratified polytrope 
whose EBK eigenfrequencies are given by equation ([H]). In order to expand our investigation to 
convective motions, we include a Doppler term in the Hamiltonian of equation (|l§|): 

u -k-v = (c 2 k 2 +uj^ (38) 

where v is the velocity field which is taken to be zero for the reference state. 

Writing the Hamiltonian in the form of equation (|35"1), we have: 



(39) 



where cq is the polytropic sound speed and the cutoff frequency of equation fl26|) has been used for 
lo c . The true definition of the cutoff frequency, given in equation (|l0|), also involves the density 
scale height H p which, in a consistent model, will be perturbed when the sound speed is perturbed. 
However, this effect should be small and we assume in this work that H p keeps its polytropic value 
of z//j,, independent of any perturbations to c. 

Although we have analytic expressions for the sound speed and velocities of the final states 
in the examples we consider below, we anticipate that in future applications we will not. In 
preparation, we discretize the propagation region and specify the sound speed and velocity 
fields at each point. The size of the discretization element was determined through numerical 
exp er iment ation . 

In order to integrate equations (|l9|), the code uses a slightly-modified version of the DE 
package of Shampine & Gordon (1975): a variable-order, variable-timestep Adams-Moulton- 
Bashforth PECE method. This method allows the user to specify tolerances for both the absolute 
and relative error. We found that the integration was sufficiently accurate and quick if both were 
set to 10" 5 . 



We have not yet addressed how to choose the switching time T of equation (37). In order to 
ensure adiabaticity, we must choose T to be much longer than the longest dynamical timescale 
associated with the system, in our case the time to complete one horizontal period. At first glance 
it may appear that, with infinite computer resources, increasingly better results may be gained 
by letting T — ► oo. This is not the case, however. Recall that the KAM theorem guarantees that 
some tori will be destroyed with any perturbation away from integrability (although most tori 
remain invariant for small perturbations). If the system switches through one of the destroyed 
tori and remains there for an extended period, the trajectory will begin to wander through phase 
space in a process called Arnold diffusion ( [Lichtenberg fc Lieberman 1983 ). The rate of Arnold 



diffusion depends on the strength of the perturbation, but it is always present. Thus, the choice of 
T is also bounded from above. Following the lead of other investigators, we choose T to be ~ 40 
horizontal periods. 
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Beginning from different initial conditions on the same invariant torus should, in principle, lead 
to identical eigenfrequencies for the final state. Our numerical experiments found a small scatter, 
presumably arising from the Arnold diffusion described above. This scatter is observed by other 
investigators and it has been found empirically, although there is some mathematical justification 
presented in Skodje & Cary (1988), that better results are obtained by implementing the procedure 
for several randomly distributed initial conditions and averaging the final eigenvalues. We typically 
average 20 such initial states in this work, a number which we find suitably constrains the result. 
A typical integration of one trajectory takes ~ 1 hour to complete on a 194 MHz SGI Power Onyx. 

Finally, we have scaled our variables in order to simplify the numerics. All lengths are quoted 
in units of 10 9 cm and all times in units of 10 2 ' 5 seconds. For example, the solar surface gravity is 
g = 2.7397 and the Sun's radius is R & = 69.599. 



4. Results 

Having outlined the method of adiabatic switching, we now apply it to several examples. 
Although we only consider fairly elementary problems in this work, adiabatic switching can in 
principle be applied to arbitrarily complex sound speed and velocity distributions. However, we 
believe that even the application to the simple systems considered here is new. The major physical 



constraint arises from the ray approximation made in §2.2 which allowed us to write the system 
in Hamiltonian form. Because of this approximation the validity of the method is restricted to 
systems where the perturbations from the reference state are large compared to the lengthscale of 
the mode, k~ x . 



4.1. Height-Dependent Sound Speed Perturbation 

We first explore an envelope with only a vertical perturbation to the sound speed and no 
velocity perturbation. The horizontal invariance of the final state implies that k x is a constant 
of the motion and hence that the system is integrable at point during the switching. Since 
all intermediate invariant tori exist, we expect adiabatic switching to work quite well. For a 
reasonably-chosen perturbation we can check our results using EBK quantization. In order to 
compare the two methods, we: 

1. Choose the initial conditions so that the system is in one of the quantized states given by 
equation (|3l|). 

2. Use the method of adiabatic switching to calculate the eigenfrequency of the new system 
using the perturbation described below. 



3. Use equation (27) to calculate numerically the eigenfrequency of the new, perturbed system 
and compare with the results from adiabatic switching. 
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Consider the simple, but not particularly realistic, case where the sound speed is close to the 
polytropic version given by equation (p6|): 

— V [l + €sin(«z)] , (40) 
lij 

where e is a small parameter characterizing the sound speed variation and k is the perturbation 
wavenumber. For the ray approximation to be valid we must have k <C k z , a requirement which 
can only be met over part of the trajectory since k z = at the upper and lower turning points of 
the cavity. This difficulty is merely the well-known breakdown of the WKB approximation near 
the classical turning points. We will only require that n <^ik z over most of the vertical extent of 
the cavity. 

The acoustic cutoff frequency of equation ( |26| ) will also change to 

a; c =(^t^)"[l + 6 S inM], (41) 



As mentioned in §3.3, we ignore any accompanying perturbation to the density scale height H 



p 



arising from the sound speed perturbation. 



We compute relative frequency shifts, 5uj/lo, via adiabatic switching and EBK quantization 
and list the results in Table |l|. Even for e = 0.16, there is no appreciable difference in the frequency 
shifts. When applying adiabatic switching to non-integrable systems, however, we expect the 
trajectory to diffuse in phase space as invariant tori disappear, and thus the computed frequency 
shifts will be less accurate as we increase the perturbation strength. 

A closer inspection of Table [l] reveals two interesting trends. First, all frequency shifts are 
positive, and second, as the final five rows demonstrate, the frequency shift is essentially linear 
in the perturbation parameter e. Following the analysis of Gough (1993), we can explore both 
of these effects by expanding the EBK quantization condition of equation (^0|). We write the 
shifted frequency as u = ujq + 8u>, where loq is the unperturbed frequency. Similarly, we have from 
equation (40) 



i , i 

2 



c = Co + 5c= \ 2 + e ' sin(Kz) ( 12) 

where cq is the polytropic sound speed. 

Applying the EBK quantization of equation ((30) to the final state yields 



22 rz% 
k z dz = I 

21 J 21 



(loq + 5uj) 2 — UJ 2 



(co + 5cf 



c h 2 

hj, r . 



2 TT 



dz = -, (43) 



where we have selected the mode with n = 1, ignored any perturbations to the cutoff frequency, 
and taken the integral over the original cavity (which is assumed to remain unperturbed). Next, 
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assuming that 5ui and 5c are small parameters, we expand and discard higher-order terms. After 
rearranging, we have 



Z-2 



Cz,0 



1 + 



r 2u2 



woo^ — <^o 1" w 

V Co 



dz 



7T 



(44) 



where /c 2j o is the positive root of equation fl29|). We expand again and cancel the zeroth-order 
term, leaving (after some algebraic simplifications) 



8u 

UJQ 



z * 5c kl 
i c k Z)0 



dz 



22 

i CqA; 2j o 



■<iz 



(45) 



Thus, in a rough sense, 5u is related to the integral of 5c over the cavity. In particular, 
since 5c as given in equation depends linearly on e, we expect the same to be true of 5u>, an 
expectation confirmed by adiabatic switching. Furthermore, since every other term in equation 
( p5| ) is non-negative, the sign of the frequency shift is determined by the sign of the integral of 5c 
over the cavity. 

Figure [2] shows the sound speed perturbation for the parameters e = 0.01 and k = 0.10. The 
dashed lines indicate the vertical extent of the unperturbed resonant cavity. The perturbation is 
positive through the cavity, leading to the upward frequency shifts displayed in Table [|. 



4.2. Cell-like Velocity Perturbation 



In both this subsection and the next we will explore simple models of a convective cell. While 



the structure of the cell in §4.2 is more realistic, its complexity hides some important physical 
results. So, as a first attempt we model a convective cell as a velocity perturbation with no 
accompanying sound speed perturbation. Our velocity field is derived from a stream function 



y x V 



/2vrx\ . /2vr(z + 6) 
e sin I — — J sin 



L 



(46) 



where L = 2ttRq is the length of the cavity, e and b are constants, and y is a unit vector 
perpendicular to the previously-defined x- and z-axes.. Although this parameterization, by 
assuming incompressibility, violates the conservation of mass flux, we are more concerned with the 
physical insight we can derive from this model than its self-consistency. Expanding equation ( [46] ) 
gives the velocity field 



■ sin 



x 

R & 



COS 



z + b 



and 



R G 



cos 







X 

R & 



sin 



z + b 

Rp) 



(47) 



Figure |3| shows the streamlines for the velocity profile given by equation ( |47| ) with 6 = 3. The 
two dashed lines near the top of the plot delineate the resonant cavity of the initial eigenstate. 
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The WKB criterion is clearly met in the vertical direction. The streamlines close at lower depths, 
which are not shown in order to better resolve the cavity. Units for the x- and z-axes are those 



given in §4.1 and the contours in each plot are linearly spaced. 



For all values of b except integer multiples of 2ttR & , the vertical velocity is unphysically 
non-zero at z = 0. However, we can see from Figure | that by varying b, the effects of propagation 
through different portions of the convective cell are easily studied. We make the offset parameter 
b non-dimensional by defining 

b + z c 



(48) 



where z c is the mean depth of the resonant cavity. The method of adiabatic switching is then 
used to investigate the dependence of eigenfrequencies on the parameters e and b. Our results are 
shown in Table 

Clearly the frequency shift depends on the strength of the perturbation in a different manner 
than the system we considered in the previous section. Instead of a linear dependence, the shift 
appears to be second-order in e. In addition, the frequency shift is always downward. We can 
explain both these effects through the mathematical framework we developed earlier. 

Our system contains no perturbation to the sound speed, but now the Doppler term of 
equation (pSl) must be included when evaluating equation (|30|): 



k z dz 



(ljq + 5u> — k • v) z 



1/2 



dz 



TT. 



(49) 



where we have explicitly kept the line integral from equation (27). Performing the same expansion 



as done in §4.1 while treating k • v as a small parameter, we arrive at 



8u) 



k- v 



rdz 



-l -1 



^0 



dz 



(50) 



Again, we find that the frequency shift depends on the integral of the perturbation over the 
cavity. Since, as can be seen from Figure ||, the ray does not interact with all of the convective 
cell, we might expect the perturbation term not to average to zero. However, in this case the 
perturbation has an almost antisymmetric effect on rays propagating with and against the flow. 
To first order these shifts cancel, leaving 5uj = 0. However, we expect a non-zero second-order 
frequency shift, as can be shown by the following thought experiment. To first order, propagating 
with and against the velocity field produces equal and opposite Doppler shifts of the ray's 
frequency. But, the ray takes more time to traverse regions where the ray travels against the flow 
(k • v < 0) than it does to cross regions where the ray and the flow are co-directional (k • v > 0). 
The relative time difference is of order v/c and the value of the frequency shift is of order —v/c, 
leading to a net downward shift of order —v 2 /c 2 . This result is not restricted to the case considered 
here, but instead applies for any perturbation where the ray spends equal times, to first order, in 
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upshift and downshift regions. In particular, sound perturbations caused by the thermal structure 
of a convective call have this effect. 

Instead of carrying out the expansion of equation (p9|) to second order in the small parameters, 



we can graphically demonstrate this result by varying the parameter <fi of equation (48). By doing 
so, the resonant cavity will overlap different regions of the convective cell, leading to different 
frequency shifts. Near the top and bottom of each cell, the flow is predominantly horizontal. 
When the resonant cavity overlaps this region the horizontal component (k x v x ) of equation (|50|) 
will produce the largest shift. Similarly, when the cavity overlaps regions of strong upflows and 
downflows, the vertical component (k z v z ) will produce the strongest influence. 

Figure |] shows the relative frequency shift plotted versus the phase (f>, for a n = 1, i = 200 
mode with z c = 0.87 and e = 0.05. The structure of the convective cells in equation (|46| ) establishes 
that the plot is periodic outside the plotted range of < <fi < ir. In this example, k x is always 
larger than \k z \ and so we expect the effects from the horizontal perturbation to be the strongest. 
The points </> = and eft = tt correspond to configurations where the resonant cavity lies near 
the vertical intersection of two cells, and hence we would expect these points to show the largest 
perturbation, which is indeed the case. The minimum frequency shift is found, as expected, when 
<j> = tt/2 since for this alignment of the cell and the resonant cavity the horizontal velocities 
sampled by the ray are smallest. 



4.3. Convective Cell 

While the cell studied in the preceding section was helpful in shaping our intuition, it was 
not particularly realistic. A more consistent model has sound speed and velocity perturbations in 
both the x and z directions. In this section we explore such a model and will show that the net 
frequency shift is still negative and second order in the strength of the perturbation. 

Our reference state remains the adiabatically stratified, plane-parallel polytrope described 
above. We take our model for the convective state from Shirer (1987). It is a simple nonlinear 
model of steady Rayleigh-Bernard convection derived from the Boussinesq equations in two 
dimensions. We work with the solution with the lowest mode numbers, two cells in the horizontal 
direction and one in the vertical. 

After converting Shirer's parameterization to the notation of this paper, we are left with three 
free parameters: the depth of the convective cell zt, the thermal diffusivity k, and its product kv 
with the kinematic viscosity. Mindful of the constraint that the lengthscales of the fluid motions 
should be much larger than the lengthscales of the oscillation, we work with the largest possible 
convection cell, one the approximate depth of the solar convection zone, zt = 0.3i? Q . 

Estimates of k and v in the solar convection zone vary widely over a number of orders of 
magnitude. Our choice for these terms, however, does not have a physical motivation. Instead, 
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we choose these constants so that the size of the perturbations is kept small. A justification for 
this assumption comes from the observation that the solar convection zone is nearly adiabatically 
stratified, hence we expect the convective motions to be perturbations on the mean sound speed 
structure. 



With these choices, the sound speed becomes 



2 9 Z , 
c = V a cos 

A* 



5111 



7T ( 1 



/3sin 



2ir ( 1 

ZT 



(51) 



where L = 2ttRq is the length of the cavity and a and (3 are functions of the product kv. The 
velocity perturbations are given by 



7 8111 



2-7TX 



cos 



vr 1 



2 



and 



6 cos 



2vrx\ 



sm 



vr 1 



(52) 



where 7 and <5 are functions of k and the product kv. The nonlinear nature of the convection arises 
from the final term in equation ([51]) which is independent of x. For a linear cell, the perturbation 
would oscillate in x. 



Specifically, 



gz T V8 a 2 + 1 1 
a = Ka — Ka r 2 



/x 7T Ra a 



(53) 



where 



2z T _ 0.3 



Ra 



7T 



g4 

KVTT 4 



and 



Ra r 



(a 2 + l) 3 



(54) 



are the aspect ratio of the convective domain, the Rayleigh number, and the critical Rayleigh 
number for the onset of convection, respectively. The parameter (3 is defined by, 



f3= 9ZT * [Ra-Ra c 



We also have 



and 



fj, it Ra 
TTy/8 



z T (a 2 + I] 

7T\/8a 
z T {a 2 + 1) 



k [Ra — Ra c ] 2 
k [Ra — Ra c ] 2 



(55) 



(56) 



(57) 



From the expressions in equations (53)— (57) we can see that a critical parameter governing 
the sizes of the perturbation is 



(Ra — Ra c ) 



g4 1 

7T 4 KV 



(a 2 + If 



(58) 



In order to keep the perturbations small, we must choose the free parameters - k and v - such 
that this expression is small, but non-zero. For instance, choosing kv = 2.27 and k = 0.1 gives 
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Ra - Ra c = 9.5 x 1(T 2 , a = 0.50, = 5.1 X 10~ 3 , 7 = 0.013 and 5 = 1.2 x 10~ 3 . Since we are 
exploring the parameter regime near the boundary of convective instability, small changes in the 
parameters can have large effects on the size of the perturbations. 

Figure || shows the streamlines for the velocity profile given by equation (52) as well as the 
isotherms for the temperature perturbation given by the last two terms of equation (|5l|). The 
standard picture of rising warmer fluid and sinking cooler fluid is seen to hold. As in Figure || the 
streamlines and isotherms close at depth and are linearly spaced. 

This convective cell can be thought of as a combination of the perturbations considered in 
the two previous subsections. Comparison of the velocity perturbations of equations (|52|) and (fl7|) 
show they are quite similar. We expect the perturbations arising from the velocity profile to be 
second-order by the argument advanced earlier. The sound speed perturbation of equation (|5l| ) is 
an interesting combination of a strong perturbation which we expect to enter the frequency shift 
at second order (the second term) due to its periodic horizontal structure and a small perturbation 
which will enter at first order (the third term). Our results show that it is difficult to separate the 
effects of these terms, at least in the parameter regime we are considering. The different runs are 
documented in Table ^ and show that the frequencies are downshifted and generally are of second 
order in the perturbation strength. 

More interesting than the frequency shifts is the fate of the raypaths under the influence 
of this perturbation. As the most realistic model of convection under consideration, we expect 
that these raypaths will most closely approximate raypaths in the solar convection zone. Recall 
from Figure [l] that the raypaths and caustics for integrable systems are quite regular. We expect 
from the KAM theorem that for small perturbations the eigenrays will remain on invariant 
tori, implying that both the caustics and the raypaths will remain similar to, although possibly 
deformed from, the integrable shape. In addition, as the strength of the perturbation increases we 
expect the tori will eventually break down and the ray will begin to explore the entire coordinate 
space. 

Figure ^ shows that both of these expectations are borne out. We have plotted sample 
raypaths for two different strengths of the convective perturbation. These raypaths were computed 
by running the raypath integration routine after adiabatic switching had completed, hence the 
frequency is constant in each panel. As in Figure |], the horizontal structure of the ray has been 
altered for legibility. The top panel shows the raypaths for a weak perturbation and confirms 
our prediction that the ray will be confined to an invariant torus. The bottom panel shows the 
result from a stronger perturbation, the path of a ray which not confined to an invariant torus. If 
this integration was continued indefinitely, the ray would eventually diffuse onto a near-vertical 
trajectory and plunge out of the bottom of the displayed propagation region. 
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5. Discussion and Future Work 

We have demonstrated how to compute frequency shifts arising from general perturbations 
with no special symmetries. From our investigations, we can draw some general conclusions 
concerning the effects of convection on p-mode eigenfrequencies. To first order, the relative 
frequency shift is roughly proportional to the integral of the perturbation over the resonant cavity 



of the eigenmode. For cases such as the sound speed perturbation considered in §4.1 where the 
perturbation is always either positive or negative, the shift is linear in the perturbation strength 
while the direction of the shift depends on its sign. 

Of more physical interest is the case where the perturbation is not monotonic, but periodic. 
Here, we have shown through physical and mathematical arguments as well as numerical methods 
that the effect of the perturbation vanishes to first order. However, a second-order perturbation 
does exist which has a pleasing interpretation of being due to the extra time spent in regions of 
counter-propagating material. It is intriguing to note that, as discussed in Gough et al. (1996), the 
structure of present solar models near the top of the convection zone yields p-mode eigenfrequencies 
which are larger than those observed by helioseismology. This results of this paper show shifts in 
the proper direction; but, we have not yet shown that the second-order effect we discuss above is 
large enough to account for these differences. 

In addition, we have demonstrated that even mild perturbations can significantly change the 
structure of an eigenmode's resonant cavity despite producing only small frequency shifts. As the 
perturbation strength increases, the ray eventually begins to move chaotically and is no longer 
confined in any obvious region. It is an open question whether this effect will persist when we 
investigate more realistic convective models. If it does, results from time-distance helioseismology 
imply that the ray parameterization does work in the Sun suggesting, perhaps, that consideration 
of average and not individual raypaths would be useful. 

We have shown that adiabatic switching is a viable method of implementing the EBK 
quantization conditions for nonintegrable systems and finding the eigenfrequencies of simple 
convective envelopes in the WKB limit. For integrable systems where the eigenfrequencies can be 
directly computed from EBK quantization, adiabatic switching agrees quite well. Furthermore, 
adiabatic switching can find eigenfrequencies of non-integrable systems which are impractical to 
obtain through direct EBK quantization. 

Adiabatic switching is an attractive method for finding eigenfrequencies, but care must be 
taken to apply it to problems for which it is valid. In the Sun, the true excited objects are global 
modes which are not necessarily well approximated by one-dimensional objects such as rays. 
Mathematically, one may approximate a ray whose dynamics are governed by equations ( |i~9| ) by 
superposing eigenmodes to form a wavepacket. But, Bogdan (1997) has shown that, due to the 
finite number of excited p modes, the smallest wavepacket which can be constructed by a coherent 
superposition in the Sun is ~ 30 Mm. This dimension, characteristic of supergranules, suggests 
a lower size limit on the structures which can be explored using adiabatic switching. Smaller 
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structures are incompletely sampled by the wavepacket, suggesting that their effects are not fully 
realized. The smallest convective structures cannot be treated at all by WKB theory and their 
effects may well be modeled by turbulent pressure ( Goldreich fc Murray 1994 , Pruzinov 1998; ). 



Provided that they are of large enough scale, convective motions of almost arbitrary 
complexity can be investigated with adiabatic switching. A number of intriguing problems remain 
in this area including the exploration of other analytic models of convective structures such as 
upwellings and thermal plumes. We are not limited to analytic models, however. Numerical 
models of convection such as recent spherical shell simulations ( |Glatzmaier Sz Toomre 1995| , Elliot 



et al. 1999 ) and convection in layers ( Brummell et al. 1998| ) can also be studied with adiabatic 



switching. 

Along with exploring convective structures, we can also investigate modifications to the 
dispersion relation of equation (|38|). A fairly simple conceptual, although perhaps computationally 
expensive, extension to this work would be to incorporate the third dimension into the plane- 
parallel geometry. The number of degrees of freedom would increase by one, however the reference 
system would remain integrable since another constant of the motion, k y , would also be added. 



Perturbations analogous to the Doppler term of equation (38) can also be included. Both rotation 



and magnetic fields can be expressed as additional terms in the dispersion relation ( Gough 1993| ) 



Power spectra of the solar oscillations show that p modes possess an intrinsic line profile, the 



source and shape of which is not well understood. Previous works ( Kumar et al. 1994 , Roxburgh 



fc Vorontsov 1995] , Rast fc Bogdan 1998 ) have treated the effects of intrinsic damping, noise, and 



source structure. Intriguingly, eigenfrequency linewidths are also generated by some of the effects 
we have explored in this paper. We ignored the inherent time-dependence of solar convection in 
this work; recall, however, that the frequency is a constant of the motion only for time-independent 
systems. Indeed, only by introducing X(t) into equation ( |35| ) could adiabatic switching work. In 
addition, the Sun is almost certainly not an integrable system, leading us to expect that a ray 
will diffuse in phase space as it propagates. Since the frequency is not a constant of the motion, 
its frequency can gradually change. We saw this effect in our work with non-integrable systems, 
but averaged results from several initial conditions to obtain a final result. Perhaps a closer 
exploration of the raw results will allow us insight into the nature of p-mode linewidths. 

Finally, as a byproduct of using adiabatic switching, we generate multiple raypaths together 
with travel times for a given envelope. This library of paths could prove useful in investigations of 
time-distance helioseismology. 
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Table 1. Eigenfrequencies for the vertical sound speed perturbation of 
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a Vertical mode number from equation (^Tj) 
b Horizontal mode number from equation ( |33|) 
c Perturbation parameters from equation ( ff0| ) 
d Relative frequency shifts from adiabatic switching 
e Relative frequency shifts from EBK quantization 

f The first three rows give unperturbed frequencies for comparison purposes 



Table 2. Eigenfrequencies for the velocity perturbation of 
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a Vertical mode number from equation ( |3~I| ) 
b Horizontal mode number from equation ( |33[ ) 
c Perturbation parameter from equation fl47| ) 
d Perturbation parameter from equation (|4^ ) 
e Relative frequency shifts from adiabatic switching 
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Table 3. Eigenfrequencies for the convective cell of §4.2 
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a Vertical mode number from equation ( |3ll) 
b Horizontal mode number from equation ( |33[ ) 
c Thermal diffusivity - see equations (|53| - 57) 
d Kinematic viscosity - see equations (p3 - 57) 



e Relative frequency shifts from adiabatic switching 
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Fig. 1. — Sample raypath in an adiabatically stratified, plane-parallel, horizontally-periodic 
envelope. The cavity has been horizontally compressed for clarity. Caustic surfaces at z = 0.32 
and z = 1.42 are visible. 

Fig. 2. — Relative perturbation to the sound speed described in equation ([40|). The dashed lines 
show the vertical extent of the ray's unperturbed cavity. 



Fig. 3. — Streamlines for the velocity perturbation of §4.2 with 6 = 3. The streamlines connect at 



lower depth and the dashed lines show the resonant cavity of the intital ray. 

Fig. 4. — Relative frequency shift versus the phase of the offset between the resonant cavity and 
the convective cell as defined in equation (| 



Fig. 5. — Upper panel: Streamlines for the convective perturbation described in §4.3. The 
streamlines connect at lower depths. The caustic surfaces of a sample raypath are shown as dashed 
lines. Lower panel: Temperature perturbation isotherms for the same convective structure. 

Fig. 6. — Upper panel: Sample raypath for a weak convective perturbation (k = 0.01, v = 227). 
As in Figure (Q), the horizontal scale has been compressed for legibility. The deformed caustics 
are clearly visible. Lower panel: Raypath for a stronger perturbation (k = 0.1, v = 22.6). The 
ray is no longer confined by caustics and will eventually explore all of coordinate space. Note the 
difference in vertical scales. 
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Relative Sound Speed Perturbation 
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Streamfunction for velocity perturbation 
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Streamfunction for convective perturbation 
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Isotherms of Temperature Perturbation 
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Raypath from Weak Perturbation 
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